Optimization of synchronization in gradient clustered networks 
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We consider complex clustered networks with a gradient structure, where sizes of the clusters are distributed 
unevenly. Such networks describe more closely actual networks in biophysical systems and in technological 
applications than previous models. Theoretical analysis predicts that the network synchronizability can be 
optimized by the strength of the gradient field but only when the gradient field points from large to small 
clusters. A remarkable finding is that, if the gradient field is sufficiently strong, synchronizability of the network 
is mainly determined by the properties of the subnetworks in the two largest clusters. These results are verified 
by numerical eigenvalue analysis and by direct simulation of synchronization dynamics on coupled-oscillator 
networks. 
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It has been recognized in biological physics that at the cel- 
lular level, information vital to the functioning of the cell is 
often processed on various networks with complex topolo- 
gies 111]. At a systems level, organizing information using 
the network idea has also become fundamental to understand- 
ing various biological functions. A key organizational feature 
in many biological systems is the clustered structure where 
biophysical and biochemical interactions occur at a hierarchy 
of levels. Examples include various protein-protein interac- 
tion networks Q2L I2D and metabolic networks |4j] . In biology 
and network science, a fundamental issue is synchronization 
llH @|. The aim of this paper is to study synchronization in 
clustered complex networks with uneven cluster-size distribu- 
tion and asymmetrical coupling. Since this type of network 
structure is also important to physical and technological sys- 
tems such as electronic-circuit networks and computer net- 
works [7, 8, 9], understanding synchronization in such net- 
works will be of broad interest. 

There has been recent effort to study synchronization in 
complex clustered networks fiol [Till . A general assumption 
in these works is that all clusters in a network are on the equal 
footing in the sense that their sizes are identical and the in- 
teractions between any pair of clusters are symmetrical. In 
realistic applications the distribution of the cluster size can 
be highly uneven. For example, in a clustered network with 
a hierarchical structure, the size of a cluster can in general 
depend on the particular hierarchy to which it belong. More 
importantly, the interactions between clusters in different hi- 
erarchies can be highly asymmetrical. For instance, the cou- 
pling from a cluster at a top hierarchy to a cluster in a lower 
hierarchy can be much stronger than the other way around. An 
asymmetrically interacting network can in general be regarded 
as the superposition of a symmetrically coupled network and a 
directed network, both being weighte d. A weighted, directed 
network is a gradient network IU2[ 11311 , a class of networks for 
which the interactions or couplings among nodes are governed 
by a gradient field. Our interest is then the synchronizability 
and the actual synchronous dynamics on complex clustered 



networks with a gradient structure. 

For a complex gradient clustered network, a key parame- 
ter is the strength of the gradient field between the clusters, 
denoted by g. A central issue is how the network synchroniz- 
ability depends on g. As g is increased, the interactions among 
various clusters in the network become more directed. From 
a dynamical-system point of view, uni-directionally coupled 
systems often possess strong synchronizability Ill4ill5ll . Thus, 
intuitively, we expect to observe enhancement of the network 
synchronizability with the increase of g. The question is 
whether there exists an optimal value of g for which the net- 
work synchronizability can be maximized. This is in fact the 
problem of optimizing synchronization in clustered gradient 
networks, and our findings suggest an affirmative answer to 
the question. In particular, we are able to obtain solid analytic 
insights into a key quantity that determines the network syn- 
chronizability. The theoretical formulas are verified by both 
numerical eigenvalue analysis and direct simulation of oscil- 
latory dynamics on the network. The existence of an optimal 
state for gradient clustered networks to achieve synchroniza- 
tion may have broad implications for evolution of biological 
networks and for practical applications such as the design of 
efficient computer networks. 

Our general setting is network with N nodes and M clus- 
ters, where n m is the size of cluster m and V m denotes the 
set of nodes it contains (m = 1, M). Each pair of nodes 
is connected with probability p s in the same cluster and with 
probability pi in different clusters, where p s > pi IToll . For 
a coupled oscillator network with arbitrary connecting topol- 
ogy, its synchronizability is determined [16] by the interplay 
between the transverse stability of the local-node dynamics 
F(x) and the eigenvalue spectrum of the coupling matrix C, 
which can be sorted conveniently as Ai = < A2 ^ • • • Sj 
An, where Ai = underlies the synchronization solution. A 
typical nonlinear oscillator in the synchronization manifold is 
transversely stable only when some generalized coupling pa- 
rameter a falls in a finite range: a G [eri, 02], which is deter- 
mined by the single-oscillator dynamics. The network is syn- 
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chronizable if all the normalized eigenvalues except Ai can be 
contained within this range: oi < eA 2 ^ • • • ^ eXjf < u%, 
where e is a specific coupling parameter. For convenience, 
we consider the following class of coupled-map networks: 
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e?-dimensional map representing the local dynamics of node i, 
e is a global coupling parameter, and H is a coupling function. 
The rows of the coupling matrix C have zero sum to guarantee 
an exact synchronized solution: x* = x^ = ... = x^ = s t . 
For certain types of oscillator dynamics and coupling func- 
tions, say, for example, the linearly coupled logistic oscilla- 
tors we are going to study in the following, <tn is sufficiently 
large 11711 . In such cases the condition eXn < &i is naturally 
satisfied and the synchronizability of network is only deter- 
mined by A2. For simplicity, we will restrict our study to such 
types of oscillator dynamics and coupling functions. 

We first develop a theory for networks consisting of two 
clusters (the theory can be generalized to multiple-cluster net- 
works). Without a gradient field, the adjacent matrix A is such 
that A 
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1 if there is a link between node i and node j, 
and Aij — otherwise. To introduce a coupling gradient 
field from cluster 1 to cluster 2, for each inter-cluster link 
i € V\ and j G V 2 , we deduce an amount g from 
A^ (corresponding to the coupling from node j to node i) 
and add it to A,; so that the total coupling strength is con- 
served. In this sense the gradient field can be said to point 
from cluster 1 to cluster 2. The coupling matrix C is defined 
as Cij — —Aij/ki, where ki = YljLi Aij is the weighted 
degree of node i, and Cu = 1. 

The eigenvalue spectra of C and of its transpose C T 
are identical. Let e 2 = (ei, e 2 , e ni , e ni+i , e N ) T 
be the normalized eigenvector associated with A 2 of C T . 



Since YljLi ^J,i = 2j=i < ^ i i = ^> e ig envectors as- 
sociated with non-zero eigenvalues of C T have zero sum: 
Ylf=i^,j — [18]. From C T e 2 = A 2 e 2 we have A 2 = 



e^C T e 2 = ^ii=i e iCij e j- F° r a clustered network, the el- 
ements in e 2 have a special distribution: « E\ for i g V\ 
and tj ~ E2 for j 6 V2 [10], where the two constant 
values Ei and _E 2 can be obtained from the normalization 
condition e 2 n e 2 = 1 and the zero-sum property. We obtain 

Ei = -\Zn 2 /(nin 2 + n\) and E 2 = \Jni/(nin 2 + nl) 
(the signs of Ei and E 2 are interchangeable since EiE 2 < 
0). This can greatly simplify the calculation of A 2 , which 

now can be written as A 2 w Z)»=i e i{(Cii + + ... + 
C ini )Ei + (C ini+ i + C ini+2 + ... + C iN )E 2 }. The non- 
zero elements in C can be calculated as follows. For i G V\, 
ki w nip s + n 2 pi(l - g), if j G Vi, Qj = -1/ki = gn, 
and there are approximately nip s non-zero elements for each 
i. If j G V 2 , we have Cy = —(1 — g)/h = gi 2 . For i G V 2 , 
h ks n 2 p s + nipi(l+g), if j G V\, Qj = -{l + g)/ki = g 2 i 
and, if j G V 2 , Cij = —1/ki = .922- Since Cu = 1, the 
calculation can be further simplified as A 2 « ^27=1 e *{-^i + 
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FIG. 1: (Color online) Theoretical contour plot of A2 in the (g,ni) 
plane, for a 2-cluster network of ni + n2 = 300 nodes. Other pa- 
rameters are pi — 0.2 and p s = 0.7. The dashed line is given by Eq. 
(O, which determines, for fixed value of m, the optimal gradient 
strength go. 
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A 2 = l + {Einigii + E 2 n 2 g 22 )p s + EiE 2 nin 2 pi(gi 2 +g 2 i). 

(1) 

In Eq. (Q~|i, the unity comes from the diagonal elements in C, 
it defines the upper limit for A 2 (this special case is associated 
with one-way coupled tree-structure networks 001). The 
second term is contributed by the intra-connection of clus- 
ter 1 and cluster 2. The last term corresponds to the inter- 
connection between the clusters. The parameter g is contained 
in these terms via g^ . For a given 2-cluster network, the opti- 
mal gradient strength go that maximizes A 2 can be determined 
by setting d\ 2 /dg = 0, which gives 



2m -N 

9o = — 77 [Ps - Pi) 

Npi 



(2) 



and niE{ + n 2 E 2 = 1 (the normalization condition), we ob- 



(Please note that in deriving go we actually get two such val- 
ues: go and g = N(p s +pi)/[(N - 2ni)pi] < -1. Since in 
our network model \g\is defined within range [0,1], the value 
g a is therefore discarded.) 

Equation (fTJ reveals some interesting features about the de- 
pendence of A 2 on key parameters of the clustered network. 
To give an example, we show in Fig. Q]a contour plot of A 2 , 
calculated using the theoretical formula Eq. ([TJ, in the param- 
eter plane spanned by ni and g, where n\ + n 2 — 300. It 
gives, for fixed value of ni, the dependence of A 2 on gradient 
strength. Since, by our construction, the gradient field points 
from cluster 1 to cluster 2, the upper half region (rti > 150) 
in Fig. Q] represents gradient clustered networks for which 
the gradient field points from the large to the small cluster. 
For any network defined in this region, for any fixed value of 
g, A 2 increases monotonically with m, indicating enhanced 
network synchronizability with the size of the large cluster. 
However, for a fixed value of n\, A 2 first increases, reaches 
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FIG. 2: (Color online) For a gradient network of two clusters with 
N = 300 nodes, numerically obtained (circles) dependence of A2 
on the strength g of the gradient field for the two cases where (a) the 
gradient field points from the larger to the small cluster (ni — 190 > 
N/2) and (b) the opposite (ni = 110 < N/2). The solid curves are 
from theory, (c) For ni = 190, actual synchronization time versus g 
for a clustered network of chaotic logistic maps. We observe a sharp 
reduction in the time as g approaches its optimal value, indicating a 
stronger synchronizability. Other parameters are pi — 0.2, p a = 0.7. 
Each point is the average of 100 random realizations. 



maximum for some optimal value of g = go, and then de- 
creases with g. The dependence of go on n\ is revealed by 
the dashed line in the figure [Eq. (0]. We see that, when the 
gradient field is set to point from the large to the smaller clus- 
ter, in order to optimize the network synchronizability, larger 
gradient strength is needed for larger difference in the cluster 
sizes. In contrast, in the lower-half of Fig. Q]where n i < n 2, 
A 2 tend to decrease as g is increased (for fixed ni) or when the 
difference between the sizes of the two clusters enlarges. This 
indicates that, when the gradient points from the smaller to 
the larger cluster, the network synchronizability continuously 
weakens as the the gradient field is strengthened. 

To provide support for our theoretical formula Eq. (Q}, we 
consider the same network in Fig. [2] and directly calculate the 
eigenvalue spectrum for a systematically varying set of val- 
ues of g. Figure [2ja) shows A2 versus g (open circles) for the 
case where the gradient field points from the large to the small 
cluster (ni = 190 > N/2) and Fig. |2|b) is for an opposite 
case (ni = 110 < N/2). The solid curves are theoretical pre- 
dictions. We observe a good agreement. To gain insight into 
the actual dynamics of synchronization on the network, we 
use the logistic map f(x) = 4x(l — x) as the local dynam- 
ics, e = 1, and choose H(x) = x as the coupling function. 
For the logistic map, we have a\ = 0.5, 02 — 1-5 11911 . We 
find numerically A at ks 1.1 < cr%. Thus the synchronization 
condition becomes A2 > <J\ = 0.5. We have calculated the 
average synchronization time T as a function of g, where T is 
the time needed to reach Y*Li \ (x l - {x)) \ / N < 5 = 10~ 5 
and (x) = YliLi x l /N (the system is considered as unsyn- 
chronizable when T > 10 4 ). As g approaches the optimal 
value go, we observe a sharp decrease in T, as shown in Fig. 
|2|c), indicating a significant enhancement of the network syn- 
chronizability. After reaching the minimum at go, the time 



increases as g is increased further, as predicted by theory. 

The theory we have developed for two-cluster networks 
can be extended to multiple-cluster networks. Consider 
a M-cluster network, where each cluster contains a ran- 
dom subnetwork. Assume the size of the clusters satisfy 
n\ > n.2 > 713 ^ • • • ^ riM, a coupling gradient 
field can be defined as for the two-cluster case. For a 
random clustered network, the weighted degree can be 
written as fe 4 ~ J2f=i A v = n ™Ps + (N - n m )p t + 
PW(T,i,n m <n, n i - E;',„ m >n,, n i') = K ™- Define g ml as 
the average value of the non-diagonal, non-zero elements 
dj. For i e V m and j € Vj, we have g mm = -1/K m , 
9ml = -(1 - g)/K m for n m > nj, g ml = -(1 + g)/K m 
for n m < m, and g ml = -l/K m for n m = raj. For 
the second eigenvector of C T , e.g. C 1 £2 — \2&2, its 
components have a clustered structure, i.e., for all i 6 V m , 
&2,i ~ E m while they may vary significantly for differ- 
ent clusters. The eigenvalue A 2 can then be expressed 

as A 2 = e 2 "C ,T e 2 = J2i,j=i e iCijej = £)i=i &i{E m + 

E m n m Ps9mrn + J2l^m E l n lPl9ml} = Em=l ?l ">^ + 

Em=i E m n m p s g mm 

+ Yji^ m E m E i n m niPig m i- Tak- 
ing into consideration the normalization condition 

eje 2 = 1, we get A 2 = 1 + £„=i E^n^p.g^ + 

Em,l=i;lj!m E m Ein m nipig m i. 

For a general multiple-clustered network, it is mathemat- 
ically difficult to obtain an analytic formula for the quantity 
E m . However, E m can be determined numerically. Once this 
is done, the general dependence of A2 on g and subsequen- 
tially the optimal gradient strength go can be obtained. In 
some particular cases, explicit formulas for E m and A2 can be 
obtained. Focusing on the role of the gradient in determining 
the synchronizability, we consider the extreme gradient case: 
9 = 1. Numerically, we find that for this case, with respect 
to the second eigenvector £2, only E\ and E% (corresponding 
to the largest and the second largest clusters) have non-zero 
values, while for all m > 2, E m — 0. From the normalization 
condition e^e.2 — 1 and the zero- sum property Ylj=i = 

(since Ylj=i = 0), we can SOIve f° r E\ and E2 as 
Ei = -\/n 2 /(nin2 + n\) and E 2 = \fn\/(nw2 + n%). 
Noticing g\2 = 0, we finally obtain 

2 2 

A 2 = 1 + E m n liPs9mrn + ^ E m E in m niPig m l 

771=1 171,1,1^171 

= 1 + (Efnlgu + E\n\g 2 2)p s + EiE 2 nin 2 pig2i-0) 

A numerical verification of Eq. (O is provided in Fig. 3(a). 
An observation is that, except for the difference in gij, Eq. 
(O has the same form as Eq. ([TJi, indicating that A2 is mainly 
determined by the first two largest clusters and it has little de- 
pendence on the details of size distributions of the remaining 
clusters. The remarkable implication is that, for different gra- 
dient clustered networks, regardless of the detailed form of the 
cluster size distribution, insofar as the two dominant clusters 
have similar properties, all networks possess nearly identical 
synchronizability. 
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FIG. 3: (Color online) (a) For a 5-cluster network (circles) and a 10- 
cluster network (squares), A2 versus m, the size of the largest cluster. 
The solid curve is from theory [Eq. J3)]. For the 5-cluster network, 
the size of the remaining clusters are 112 = 200, 723 = 50, «4 = 30, 
ng = 20. For the 10-cluster network, we have 71,2 = 200, 113 to n\o 
are 90, 80, 70, 60, 50, 40, 30, 20, respectively. Other parameters 
are pi — 0.15 and p 3 — 0.7. For n\ < n-2, the gradient is actually 
from cluster 2 to cluster 1. Each point is the average result of 100 
network realizations, (b) For a "cortico-cortical network" of the cat 
brain, numerical results of the dependence of A2 on gradient strength 
g. Synchronization is optimized for go ~ 0.55 



The model of gradient clustered network we have investi- 
gated here is different to the asymmetrical network models in 
literature. In Ref. ilO, EH], asymmetrical couplings have 
been employed to improve network synchronization and it is 
found that, for non-clustered networks, synchronization is op- 
timized when all nodes are one-way coupled and the network 
has a tree-structure 11411 . Different to this, in our model asym- 
metrical couplings are only introduced to inter-cluster links, 
while couplings on intra-cluster links are still symmetrical. 
This special coupling scheme induces some new properties 
to the functions of the gradient. Firstly, increase of gradient 
will not monotonically enhance synchronization. That is, di- 
rected coupling between clusters, i.e. g = 1, is not always 
the best choice for synchronization. In many cases the opti- 
mal gradient stregth go is some value between and 1, while 
the exact value is determined by the other network parameters 
[Eqs. ( [113) 1. Secondly, the direction of gradient can not be 
arranged randomly, it should be always pointing from large to 
small clusters. Finally, in the case of g — 1, network syn- 
chronizability is still related to the network topology, i.e. by 
the topology of the first two largest clusters; while for non- 
clustered network, synchronizability is only determined by the 
local dynamics fl4ll . 

Can synchronization optimization be expected in realistic 
networks? To address this question, we have tested the syn- 
chronizability of a "cortico-cortical network" of cat brain, 
which comprises 53 cortex areas and about 830 fiber con- 
nections of different axon densities l20ll . The random and 
small-world properties of this network, as well as its hierar- 
chical structure, have been established in several previous pa- 
pers izih . According to their functions, the cortex areas are 
grouped into 4 divisions of variant size: 16 areas in the vi- 
sual division, 7 areas in the auditory division, 16 areas in the 



somato-motor division, and 14 areas in the frontolimbic di- 
vision. Also, by the order of size, these divisions are hierar- 
chically organized lEoll . With the same gradient strategy as 
for the theoretical model, we plot in Fig. 3(b) the variation of 
A2 as a function of the gradient strength. Synchronization is 
optimized at gradient strength about g a w 0.55. An interest- 
ing finding is that the actual average gradient of the real net- 
work, g ave ss 0.37 12211 . is deviating from the optimal gradient 
g , indicating a strong but non-optimized synchronization in 
healthy cat brain. 

While our theory predicts the existence of a gradient field 
for optimizing the synchronizability of a complex clustered 
network, we emphasize that the actual value of the optimal 
gradient field may or may not be achieved for realistic net- 
worked systems. Due to the sophisticated procedure involved 
to determine the optimal gradient strength and the actual value 
for a given network, their numerical values can contain sub- 
stantial uncertainties. A reasonable test should involve a large 
scale comparison across many networks of relatively similar 
type (say, many different animals), hopefully demonstrating 
some kind of correlation between the optimum gradient and 
the observed values. Furthermore, such a test would include 
a sense of how large the difference between the optimum and 
observed is. Due to the current unavailability of any reason- 
able number of realistic complex, gradient, and clustered net- 
works, it is not feasible to conduct a systematic test of our 
theory. (As a matter of fact, we are able to find only one real- 
world example of gradient clustered network, the cat-brain 
network that we have utilized here.) It is our hope that, as 
network science develops and more realistic network exam- 
ples are available, our theory and its actual relevance can be 
tested on a more solid ground. 

In summary, we have uncovered a phenomenon in the syn- 
chronization of gradient clustered networks with uneven dis- 
tribution of cluster sizes: the network synchronizability can be 
enhanced by strengthening the gradient field, but the enhance- 
ment can be achieved only when the gradient field points from 
large to small clusters. We have obtained a full analytic theory 
for gradient networks with two clusters, and have extended the 
theory to networks with arbitrary number of clusters in some 
special but meaningful cases. For a multiple-cluster network, 
a remarkable phenomenon is that, if the gradient field is suffi- 
ciently strong, the network synchronizability is determined by 
the largest two clusters, regardless of details such as the actual 
number of clusters in the network. These results can provide 
insights into biological systems in terms of their organization 
and dynamics, where complex clustered networks arise at both 
the cellular and systems levels. Our findings can also be use- 
ful for optimizing the performance of technological networks 
such as large-scale computer networks for parallel processing. 
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